/************************************************************************
*************************************************************************
*************************************************************************

simulate_nonlinearsim_rev.do

Runs 5 simulations meant to push the autoregressive estimator to the 
breaking point. 

 
*************************************************************************
*************************************************************************
************************************************************************/

args jobid sdset

di "ARGS: `jobid'  | `sdset'"

local logon 1

clear
clear matrix
set more off


if mi("`sdset'") {
	local sdset 1
}


//Store current state of the random number generator
global curseed = c(seed)

//Load random seeds
do "do/rngseeds`sdset'.do"

//Load globals
do "constants/simglobals.do"

if $ncores >  1 {
	set seed ${seed`jobid'}
}

	
//Set maximum number of iterations
global prevmaxiter = c(maxiter)
set maxiter 1000


if `logon' == 1 {
	cap log close

	cap mkdir "logs"
	local tmp = string(date("`=c(current_date)'", "D M Y"),"%tdCCYYNNDD")
	cap mkdir "logs/`tmp'"
	local tmp2 = clock("`=c(current_time)'","hms")
	log using "logs/`tmp'/`tmp2'simulate_nonlinearsim_rev`jobid'"
}



//Targets for fraction of firms constrained
global consttargets 60 20 1


global dobs 	1

//Do we want to try the ar_select methods?
//global arselect ar_dem_select


//Name of folder where we'll save the results
global savename sweeps_nlsim_rev

/************************************************************************
*************************************************************************
Part 1: Helper Programs
************************************************************************
************************************************************************/

//Load data generator
do "do/programs/gendatsyn.do"

//Load GNR program
do "do/programs/gnr_jpe.do"

//Load additional testing programs
do "do/programs/structural_idtests.do"

//Load ar estimation programs
do "do/programs/ar_functions.do"

//Load sharetest program
do "do/programs/sharetest.do"

/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end







/************************************************************************
Program: truemoments
Purpose: Computes true moments (must tell it the true production function)

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop truemoments
program truemoments
	syntax name(name=name), pm(real) pk(real) pl(real) [  mu(real 0) sc(real 0) ]
	
	/************************
	Compute marginal products
	***********************/
	
	/****** CES ********/
	if "`name'" == "ces" {
		//For CES they need to pass mu
		
		if mi("`mu'") {
			error "Must pass mu (elasticity of substitution) for a CES production function"
		}
		
		if mi("`sc'") {
			local sc = 1-`pm'
		}
		
		/* There seems to be an error here
		//Marginal product of K and L
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x''*(1-`pm') * Y^(  (1-`pm'-`mu')/(1-`pm')  ) * exp(omega + eps)^( `mu'/(1-`pm') ) / `X'^( 1 - `mu' )
		}
		*/
		
		tempvar I
		gen `I' = `pk'* K^`mu' +   `pl' * L^`mu'
		
		foreach X of varlist K L {
			local x = lower("`X'")
			
			
			
			//This is actually the elasticity times Y/X
			gen MP_`X' =  `sc' *`p`x'' * `X'^`mu' / `I'  * Y/`X'
		}
		
		//Marginal product of M is just the cobb-douglas expression
		gen MP_M = `pm' * Y/M
		
		di "This is CES..."
		
		
	}
	
	/****** Quadratic ********/
	else if "`name'" == "quad" {
		//Marginal product for K and L is just param * Y/X
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}
		
		//Marginal product of M is slightly different
		gen MP_M = exp(omega + eps) * K^`pk' * L^`pl' * (`pm' - 2*M)
	}
	
	/****** Cobb-Douglas ********/
	
	else {
		//Marginal product is just param * Y/X
		foreach X of varlist K L M {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}	
	}
	
	
	
	
	/************************
	Moments: Elasticities and Sdev of marginal products
	***********************/
	
	//Elasticity = MP_X * X/Y
	foreach X of varlist K L M {
		local x = lower("`X'")
		
		//This may vary by firm
		gen elas`x'_tru 	= MP_`X' * `X'/Y 
		
		//Exclude the top and bottom 1 percent
		bysort t: cumul MP_`X', gen(MPpct)
		
		bysort t: egen sdevMP`x'_tru 	= sd(MP_`X') if MPpct > .01 & MPpct < .99
		
		//Clean up
		drop MP_`X' MPpct
	}

	
end






/************************************************************************
Program: getconstraints
Purpose: Sets a global equal to the proper values of the constraint set

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraints
program getconstraints
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist -1(.2)20 {
			spec`spec' 2^`const'

			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end


/************************************************************************
Program: getconstraintsnl
Purpose: Sets a global equal to the proper values of the constraint set
		This version uses the nl specifications

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraintsnl
program getconstraintsnl
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist -1(.2)20 {
			specnl `spec' 2^`const'

			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end


/************************************************************************
Program: gnr_wrap
Purpose: A wrapper for the gnr program (for use in bootstrapping)
************************************************************************/
cap program drop gnr_wrap
program gnr_wrap
	args spec laginst
	
	cap confirm var logtmp_y
	if _rc == 0 {
		//drop y k l m
	}
	else {
		rename (y k l m) logtmp_=
	}
	
	gen share = pricem*M/Y
	//Assume we know when the true process is linear or nonlinear
	if abs(`spec'-6)<.01 {
		//cap gnr y k l m "GMM" "keep" pricem nonlin
		
		if !mi("`laginst'") {
			cap gnr_jpe Y L K M share, `laginst'
		}
		else {
			cap gnr_jpe Y L K M share
		}
	}
	else {
		cap gnr_jpe Y L K M share, linear `laginst'
	} 
	
	if _rc != 0 {
	
		cap gen elasm_gnr 	= .
		cap gen elask_gnr 	= .
		cap gen elasl_gnr 	= .
		cap gen ar_gnr 		= .
	
	}

	drop share

end




/************************************************************************
Program: makebssamp
Purpose: Makes a single bootstrap dataset
************************************************************************/

cap program drop makebssamp
program makebssamp

	/******
	Step 1: Prepare the original dataset for merging
	******/
	
	//Make a household identifier that is numbered consecutively
	egen idcon = group(id)
	
	//Get the max and min value of IDs
	sum idcon
	local idmax = r(max)
	local idmin = r(min)
	
	//Get the max and min values for time
	sum t
	local tmin = r(min)
	local tmax = r(max)
	
	//Number of observations
	egen tag = tag(idcon)
	count if tag
	local nobs = r(N)
	drop tag
	
	
	//Save the source
	tempfile bssource
	save `bssource'
	
	
	/******
	Step 2: Construct framework of BS sample
	******/
	
	//Draw sample
	clear
	set obs `nobs'
	gen idcon = `idmin' + int( runiform()*(`idmax'-`idmin' + 1) )
	
	//The bootstrap household identifier
	gen id_bs = _n
	
	//Add ts
	expand `=`tmax'-`tmin'+1'
	bysort id_bs: gen t = _n + `tmin'-1
	
	tab t
	
	/******
	Step 3: Bring in data from source
	******/
	merge m:1 idcon t using `bssource'
	
	//These are households in the source that were not sampled
	drop if _m == 2
	
	//These are ts in the for which a sampled household is unobserved
	drop if _m == 1
	drop _m
	
	xtset id_bs t

end






/************************************************************************
Program: bspctiles
Purpose: Bootstraps the percentiles of the rk-statistic (and the vcov of the Lambda stat)
************************************************************************/

cap program drop bspctiles
program bspctiles
	//Number of reps, stub of the estimator (gnr or ar_dem), specification,
	//and a specifier for whether we should save or not
	args nreps spec savepath
	
	/******
	Step 0: Estimate true parameters
	******/
	preserve
	sort id t
	ar_dem
	local rhoest = _b[/rho]
	foreach inst of varlist k-mXm {
		gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
		gen INST_`inst' = l.`inst'
	}
	
	gen yrho2 = y - `rhoest'*l.y
	
	foreach endvar of varlist ENDO2_* {
		reg `endvar' INST_* l2.m l2.k l2.kXm
		estimates store `endvar'
	}
	
	ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id) wald fullrank
	local nel = rowsof(r(lab))
	
	local listel
	forvalues i=1/`nel' {
		local listel `listel' lab`i'
	}
	
	
	//ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id)
	restore
	//gen widstat_iv = e(widstat)

	/******
	Step 1: Prepare a dataset for the bootstrapped estimates
	******/
	cap postclose bsests
	
	tempfile bsests
	postfile bsests widnull `listel' using `bsests'
	
	
	/******
	Step 2: Prepare the original dataset for bootstrapping
	******/
	tempfile original
	save `original'
	
	
	
	
	/******
	Step 3: Bootstrap
	******/
	forvalues b=1/`nreps' {
		preserve
		makebssamp
		drop id idcon
		sort id_bs t
		
		//Preliminary construction
		ar_dem
		local rhoest = _b[/rho]
		foreach inst of varlist k-mXm {
			gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
			gen INST_`inst' = l.`inst'
		}
		
		//Rank test, etc.
		ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id_bs) wald fullrank
		matrix l = r(lab)
		local passlist
		forvalues i=1/`nel' {
			local passlist `passlist' ( `=el("l",`i',1)' )
		}
		
		//Nullify
		foreach endvar of varlist ENDO2_* {
			estimates restore `endvar'
			predict tmp
			replace `endvar' = `endvar' - tmp
			drop tmp
		}
		
		gen yrho2 = y - `rhoest'*l.y
		ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id_bs)
		local widnull = e(widstat)
		

		post bsests (`widnull') `passlist'
		restore
	}
	
	postclose bsests
	
	/******
	Step 4: Compute standard errors
	******/
	
	use `bsests', clear
	
	if !mi("`savepath'") {
		save "`savepath'", replace
	}
	
	local count 0
	foreach p of numlist 10 5 1 .5 {
		local ++count
		egen tmp = pctile(widnull), p(`=100-`p'')
		sum tmp
		local p`count' = r(mean)
		drop tmp
	}
	
	cap matrix drop labv
	corr lab*, cov
	matrix labv = r(C)
	

	
	use `original', clear
	
	foreach p of numlist 1/`count' {
		gen widstatp_`p' = `p`p''
	}
	
	matrix veclabv = vec(labv)'
	svmat veclabv


end



/************************************************************************
*************************************************************************
Part 2: Scenarios
************************************************************************
************************************************************************/


/************************************************************************
Fixed parameters
************************************************************************/

//List of constraints (in log base 2)
//local conlist -3(2)7

//Parameters of the production function
global pk 		= .11
global pl 		= .28
global pm 		= .67


//global quad 	= 10


//CES, almost cobb-douglas
global mu 		= .2   //Implies an elasticity of substitution equal to 1.25
global cespk 	= .064
global cespl 	= .33
global sc 		= .39

//CES, not very cobb-douglas
global mu2		= .8   //Implies an elasticity of substitution equal to 5
global cespk2	= .015
global cespl2	= .95
global sc2		= .39

////Standard deviations

//Anticipated shock
global sigeta1 	= .0925324 
global sigeta2 	= ${sigeta1} * 2

//Unanticipated shock
//global sigeps 	= 1.07
global sigeps 	= .2542193



//Number of firms
global nfirmnorm 	2613
global nfirmsmall 	1000


//Number of years
global nyearnorm 	7
global nyearshort 	4

//Persistance of productivity
global rhocal .7709392

global nl1 ${rhocal}
global nl2 1.1
global nl3 1



/*********************
Production Functions
**********************/

/*
The DGP program can handle production functions of the form

X(K,L)*M^pm

OR

X(K,L)*(pm1*M - M^2)

where X(K,L) is arbitrary. This works because in both cases optimal M depends on
only X(K,L), but not K and L separately.
*/

global XKL1 "gen X_kl = K^${pk} * L^${pl}"
global XKL2 "gen X_kl = (  ${cespk}*K^${mu} + ${cespl}*L^${mu} ) ^ ( ${sc}/${mu} )"
global XKL3 "gen X_kl = (  ${cespk2}*K^${mu2} + ${cespl2}*L^${mu2} ) ^ ( ${sc2}/${mu2} )"


/*********************
Productivity processes
**********************/
global nlomproc "86.30933 - 7.081727 - 33.35803*(l.omega + 7.081727) + 4.568181*(l.omega + 7.081727)^2 - .2030311*(l.omega + 7.081727)^3 "



/************************************************************************
Specifications
************************************************************************/
 
//Arguments to the DGP: klfunc mfunc constraintmult rho N T pm sigeps sigeta pricemsd


cap program drop spec1
//Baseline
program spec1
	args const 
	

	
	gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) 
	
	truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end



cap program drop specnl

//Nonlinear transition function
program specnl
	args n const

	gendatsyn  =`const' , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) ///
							nlprod("${nl1}*l.omega - `n'*${nl2}*l.omega^2 + `n'*${nl3}*l.omega^3")
	
	truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end



/************************************************************************
*************************************************************************
Part 3: Sweep
************************************************************************
************************************************************************/
local nreps = ${B} //ceil(${B}/${ncores})

cap mkdir "dta/simulations"
cap mkdir "dta/simulations/${savename}"

if ${dobs} == 1 {
	cap mkdir "dta/simulations/${savename}_bs"
}

if ${ncores} > 1 {
	cap mkdir "dta/simulations/${savename}/`jobid'_`sdset'"
	cap mkdir "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	local savefol "dta/simulations/${savename}/`jobid'_`sdset'"
	local bsfol "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	di "*******************"
	di "I am job `jobid'"
	di "*******************" _n _n _n
}
else {
	local savefol "dta/simulations/${savename}"
	local bsfol "dta/simulations/${savename}_bs"
}



/*
Estimate time to complete
How many operations?
*/
local speclist 0 2 4 
local nspecs = wordcount("`speclist'")

local totops = `nspecs' * wordcount("${consttargets}") * `nreps'
local opscom = 0
local esttime "???"



//Loop through specifications
foreach spec of numlist `speclist' {

	cap mkdir "`savefol'/spec`spec'"
	if ${dobs} == 1 {
		cap mkdir "`bsfol'/spec`spec'"
		local bshere "`bsfol'/spec`spec'"
	}
	
	di _n "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"

	
	//Get the constraints for this specification
	getconstraintsnl `spec'
	local conlist "${const`spec'}"
	local ind = 0
	
	foreach const of numlist `conlist' {
		di "-----------------------------------------------------------------------------------------"
		quietly: local ++ind
		quietly: local fc : word `ind' of $consttargets
		
		forvalues b=1/`nreps' {
			//Version of the constraint to display
			di "On Specification `spec'" _col(25) "Constraint = `const'" _col(70) "Rep = `b'" _col(85) "Est. Time: `esttime' min"
		
			quietly {
				timer on 1
		
				/********
				Generate data and compute true moments
				********/
				specnl `spec' 2^`const' 
				order elas* sdevMP*
				
				
				
		
				
				/*************************************
				Step 1: Apply Methods
				*************************************/
				
				
				//Construct second-order terms
				order k l m, last
				product k l m
				
				
				/*********
				Share test
				*********/
				gen logshare = log( pricem*M/Y )
				sort id t
				//sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(4)
				sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(2) lagdeg(1)
				
				local stub idtest_share
				gen `stub'_pr2	= r(r2)
				gen `stub'_F	= r(F)
				gen `stub'_pval	= r(pval)
				
				
				
				sharetest logshare k l m, test( l2.( c.k##c.m ) ) mmerr deg(2)
				
				local stub idtest_mmshare
				gen `stub'_pr2	= r(r2)
				gen `stub'_F	= r(F)
				gen `stub'_pval	= r(pval)
				
				
				/********
				WIDstat test
				********/
				widstat_iv_old
				widstat_iv

				
				/********
				Get estimates
				********/
				
				foreach estimator in ar_dem ${arselect} gnr {
					sort id t
					`estimator'_wrap `spec' 

				}
				
					
				
				
				//Compute standard deviation of marginal products
				//foreach method in ar_std ar_dem gnr {
				foreach method in ar_dem gnr {
					foreach var of varlist K L M {
						local sm = lower("`var'")
						gen MP_`sm' = elas`sm'_`method' * Y/`var'
						
						//Exclude the top and bottom 1 percent
						bysort t: cumul MP_`sm', gen(MPpct)
						
						bysort t: egen sdevMP`sm'_`method' = sd(MP_`sm') if MPpct > .01 & MPpct < .99
						drop MP_`sm' MPpct
					}
				}
				
				
				sort id t
				cap drop k l
				rename logtmp_* *
				order k l m kXk kXl kXm lXl lXm mXm, last

				
				/********
				Form observation
				********/
				if ${dobs} == 1 {
					bspctiles 100 1 "`bshere'/`fc'_`b'.dta"
				}
		
				
				
				keep elas* ar_* sdevMP* constrained idtest*  dofadj widstat_iv* widstatp_* labest* veclabv*
				
				
				collapse _all
				
				tempfile `b'
				save ``b''
				
				
				/********
				Estimate time remaining
				********/
				timer off 1
				timer list
				local opscom 	= `opscom'+1
				local timesofar = r(t1)
				local avgtime	= `timesofar'/`opscom'
				local timerem 	= `avgtime' * (`totops'-`opscom') / 60
				local esttime 	= string(`timerem', "%9.2f" )
				
			}
			
		
		}
		
		quietly {
			use `1'
			forvalues k=2/`nreps' {
				append using ``k''
			}
			
			
			
			//if there's a negative, change it to an underscore
			gen const 	= `const'
			gen fraccon	= `fc'
			gen spec = `spec'
			

			save "`savefol'/spec`spec'/const`fc'.dta", replace
			
			
		}
	}
	
}


//Restore default max number of iterations for GMM
set maxiter ${prevmaxiter}


if `logon' == 1 {
	cap log close
}
